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ABSTRACT 

We combine fast radiative transfer calculations with high resolution 
hydrodynamical simulations to study an epoch of early hydrogen reionisation 
by primordial stellar sources at redshifts ^5 <^ z ^ 30. We consider the 
implications of various local and global feedback mechanisms using a set of 
models which bracket the severity of these effects to determine, qualitatively, 
how they may have influenced the global star formation rate and the details of 
hydrogen reionisation. With relatively conservative assumptions, most of our 
models suggest that population III star formation proceeds in a self-regulated 
manner both locally and globally and, for a conventional ACDM cosmology, can 
significantly reionise the intergalactic medium between 15 ^ z ^ 20 as long as a 
large fraction of ionising photons can escape from these earliest galaxies. 

We then combine these results with our earlier work focusing on the role of 
population II stars in galaxies with virial temperatures ^ 10^ K at redshifts 
5 ^ 2; 20. Hence, we construct a complete reionisation history of the Universe 
which matches the Thomson optical depths as measured by the WMAP 
satellite as well as the evolution of the Gunn Peterson optical depth as seen 
in the absorption spectra of the higest redshift quasars. We find that even 
with conservative estimates for the impact of negative feedback mechanisms, 
primordial stellar sources contribute significantly to early reionisation. Future 
observations of a Thomson optical depth of Xg ^ 0.13 would bolster the claim 
for the existence of population III stars similar to the ones studied here. 
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1. INTRODUCTION 

The reionisation history of the intergalactic medium (IGM) holds many important 
clues about the onset of structure formation and the nature of the first luminous sources. 
Advances made in the past few years are making it possible to directly probe this history. 
Spectroscopic observations of distant quasars by Fan et al. (2000) and Djorgovski et al. 
(2001) revealed that the IGM was highly ionised at 2; ~ 6. Subsequent high resolution 
studies of three 2; > 5.8 quasars by Becker et al. (2001) showed a significant increase in 
Lya absorption from redshift 5.5 to 6.0 with the first possible detection of a Gunn- Peterson 
trough in the spectrum of the z = 6.28 quasar SDSS 1030.10+0524. More recently, White 
et al. (2003) presented an additional high resolution spectrum of a quasar at z = 6.37 (Fan 
et al. 2003) which also appears to have an extended Gunn-Peterson trough. Together, 
these results indicate that the Universe was reionised only a short time earlier, implying a 
relatively late epoch for complete reionisation. 

Results from the Wilkinson Microwave Anisotropy Probe (WMAP) satellite (Kogut et. 
al 2003; Spergel et al. 2003) have challenged this scenario by suggesting that much of the 
IGM was ionised earlier than inferred from the SDSS quasars. In particular, the measured 
correlation between polarisation and temperature yields an electron optical depth to the 
CMB surface of last scattering of Tg = 0.17, corresponding to an instantaneous reionisation 
epoch Zr — 17. The uncertainty associated with this measurement depends on fitting all 
parameters concerned with the TT power spectrum and the TE cross power spectrum. 
After a careful analysis, Kogut et al. (2003) have reported a 68% confidence range of 
0.13 < Tg < 0.21 corresponding to an instantaneous reionisation epoch 14 < Zj. < 20. 
In contrast, if reionisation was completed around 2; ~ 8 — 10, as is suggested by the 
SDSS quasars, the electron optical depth would be only Tg ~ 0.05 — 0.06. Although the 
discrepancy is only at the ~ 3a level, it is timely to explore the connection between the two 
sets of observations. 

In a previous paper (Sokasian et al. 2003; hereafter Paper I), we used a high resolution 
cosmological simulation in conjunction with a fast 3D radiative transfer code to study how 
stellar sources similar to those seen in local galaxies, i.e. population II type, contributed 
to the reionisation history of the IGM in the redshift interval 6 ^ 2; ^ 20. This work was 
similar to other related numerical studies (e.g., Gnedin 2000; Ciardi et al. 2000; Razoumov 
et al. 2002; Ciardi et al. 2003) except that source intensities were derived directly from 
intrinsic star formation rates computed in the underlying hydrodynamical simulations. In 
addition, the high mass resolution of the simulation allowed us to include sources down to 
M ~ 4.0 X 10^ Mq in a 10.0 h^^ Mpc comoving box, which proved to have a significant 
impact on the reionisation process. Our results showed that the star-formation history 
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inferred by Springel & Hernquist (2003) based on detailed hydrodynamic simulations and 
estimated analytically by Hernquist & Springel (2003) led to a relatively late reionisation 
epoch near 2; ~ 7 — 8 for a plausible range of escape fractions. Additionally, we generated 
synthetic spectra from a range of models with different escape fractions and showed that 
values in the range 10-20% led to good statistical agreement with observational constraints 
on the neutral fraction of hydrogen at 2; ~ 6 derived from the z = 6.28 SDSS quasar of 
Becker et al. (2001). However, the relatively late reionisation epoch predicted by these 
models cannot account for the high electron scattering optical depth inferred from the 
WMAP measurements. 

In an effort to study how this difference could be reconciled, we employed heuristic 
models with an evolving ionising boosting factor and found that in order to simultaneously 
match the WMAP and SDSS constraints, one requires an /esc = 0.20 model with a boosting 
factor that rises from unity at 2; ~ 6 to ;^ 50 near z ~ 18. In this picture, large boosting 
factors mean that the stellar production rate of ionising photons is much higher than that 
which would be produced according to our model assumptions. One mechanism that would 
boost the production rate of ionising photons would be if the stellar IMF evolved so that 
it became increasingly top-heavy at high redshift. While possible, this scenario would still 
result in an IGM that is highly ionised ai z < 13, a condition which may be inconsistent 
with the thermal history of the IGM. More specifically, Hui & Haiman (2003; see also 
Theuns et al. 2002) have shown that as long as the Universe is reionised before z — 10 and 
remains ionised thereafter, the IGM would reach an asymptotic thermal state which is too 
cold compared to the temperature inferred from Lya forest observations at 2; ~ 2 — 4. This 
conclusion, however, relies on assumptions regarding the ionising spectrum. For example, 
if the ionising spectrum is harder than the typical quasar spectrum above the H I and He 
I ionisation thresholds, then order unity changes in the ionisation fraction of H I and He I 
are necessary in the range 6 < z < 10 in order to sufficiently heat up the IGM within the 
observation constraint at 2; ~ 4. However, if the spectrum is relatively soft (more typical 
of population II stars) then this condition is avoided as long as He II reionisation occurs 
relatively late near z ~ 4 (see, for example, Sokasian et al. 2002). 

Another difficulty with boosting the production rate of ionising photons, as emphasised 
by e.g. Yoshida et al. (2003a,e) is that the IMF would need to remain top heavy for many 
generations of high mass stars. These stars have typical hfetimes ~ a few million years. 
However, the total elapsed time between 2; = 18 and 2; = 10 is more than a hundred times 
longer; therefore, it would be necessary for massive star formation to proceed efficiently 
for more than 100 generations. Studies by e.g. Omukai (2000), Bromm et al. (2001), and 
Schneider et al. (2002) indicate that very massive stars can form only if the metallicity of 
the gas lies below a critical threshold, Z ^ 10~^^^ Zq. It appears problematic for the gas in 
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star-forming halos to remain sufficiently chemically pristine for massive star formation to 
continue for over 100 stellar generations. 

The combined observational constraints summarised above suggest that the reionisation 
history of the Universe was complex, perhaps having had multiple reionisation epochs. 
Such a scenario has been considered by Cen (2003b) and Wyithe & Loeb (2003b) who 
pointed out the possibility that the Universe was first reionised by an early generation 
of massive, metal-free (population III) stars, but that then much of the IGM recombined 
once these stars were no longer able to form, and the Universe was subsequently reionised 
a second time by the next generation of population 11 stars. In hierarchical models of 
structure formation, these "first (metal-free) stars'' are thought to form through molecular 
hydrogen (H2) cooling in halos with masses larger than ~ 7 x 10^ Mq at redshifts as high 
as 2; ~ 20 — 30, provided that the halos are dynamically cold and contain sufficiently larger 
numbers of hydrogen molecules (Yoshida et al. 2003c) . 

Recent progress in simulating the formation of the very first objects in the universe 
(Abel et al. 1998; Bromm, Coppi & Larson 1999) has shown that the first luminous objects 
to form in CDM models are stars. In addition, because the first H II regions are capable of 
driving out most of the gas from the halo, it is unlikely that massive seed black holes are 
able to accrete within a Hubble time after they formed (Whalen, Abel, & Norman 2003). 
The notion that these population III stars played a role in cosmological reionisation in 
Cold Dark Matter (CDM) models has been discussed widely (e.g., Couchman & Rees 1986; 
Barkana & Loeb 2001; Mackey, Bromm, & Hcrnquist 2003; Venkatesan, Tumlinson & ShuU 
2003; Cen 2003a; Wyithe & Loeb 2003a; Yoshida et al. 2003a,b). Simulations indicate 
that these stars which form in "mini-halos" should be massive ^ 100 Mq (Abel, Bryan, & 
Norman 2000, 2002; Bromm, Coppi & Larson 2002). Whereas Abel et al. (2002) predict 
one massive star per halo Bromm et al. (2002) suggest that multiple stars can be produced 
in each dark matter halo. We discuss both scenarios in this paper. Such population III 
stars are efficient UV emitters because of their approximately constant effective surface 
temperatures near 10^ K with lifetimes of ~ 3 x 10^ yr, independent of their mass (e.g. 
Bromm, Kudritzki & Loeb 2001; Schacrcr 2002 for recent discussions). In fact, these stars 
arc capable of producing up to ~ 30 times more ionising photons per baryon compared to 
population II stars with a Salpeter IMF. 

Although the impressive yield of ionising photons from population III stars appears 
promising in terms of a possible solution to the early reionisation epoch suggested by 
WMAP, there are still many questions remaining as to whether enough of these first 
generation stars could have formed to cause cosmic reionisation. In general, semi-analytic 
modehng has been used to address these questions (see, e.g., Wyithe & Loeb 2003a,b; 



- 5 - 



Cen 2003a,b); however, the predictions of these models are often uncertain because of 
the relatively crude assumptions used to relate luminous objects to the dark matter halos 
which host them. For example, Yoshida et al. (2003c) have shown that simple criteria 
based on only the mass or virial temperatures of halos can fail by an order of magnitude 
or more because of dynamical heating by mergers. Additionally, the effects of feedback 
processes are very difficult to assess in semi-analytic models. The UV radiation in the 
Lyman-Werner (LW) band (11.26-13.6 eV) from these massive stars can photo-dissociate 
the fragile H2 molecules which are responsible for their formation. This negative feedback 
can significantly limit population III star formation and thereby substantially reduce the 
effective ionising fiux in the early Universe. Conversely, it has been argued that feedback 
from X-rays can promote H2 production by boosting the free electron fraction in distant 
regions (Haiman, Rees & Loeb 1996; Oh 2001). Efi^orts to explore feedback cfi^ccts on 
cooling and collapsing gas have been pursued by many groups (e.g., Dckel & Rees 1987; 
Haiman, Rees, & Loeb 1997; Haiman, Abel, & Rees 2000). However, these studies have 
primarily relied on semi-analytic methods which assume spherically symmetry and ignore 
details regarding formation dynamics and locations of halos in complex density fields. 

A more realistic analysis of the problem requires detailed numerical simulations capable 
of capturing the physics of early structure formation and related feedback effects. This 
was first attempted by Machacek, Bryan, & Abel (2001; see also Machacek, Bryan, & 
Abel 2003) who employed an Eulerian adaptive mesh refinement simulation to investigate 
the quantitative effects of a background radiation field on the subsequent cooling and 
collapse of high-redshift pre-galactic clouds. Ricotti et al. (2002a,b) followed up on the 
numerical approach by including radiative transfer effects and found that the formation 
of small halo objects is not significantly suppressed by the dissociating background. More 
recently, Yoshida et al. (2003c) used high resolution cosmological simulations which 
included a careful treatment of the physics of chemically pristine gas and the impact of 
a dissociating UV background to study the statistical properties of primordial halos in 
cosmological volumes. In a subsequent study, Yoshida et al. (2003b) coupled radiative 
transfer calculations to their analysis to investigate reionisation by an early generation of 
stars. More specifically, they identified plausible sites of star formation in their cosmological 
volume where they then placed massive stars with M — 300 M0 under the usual "one star 
per halo" assumption. They ran the multi-source radiative transfer code described in Paper 
I on a 200^ grid with the assumption that source lifetimes are 3 million years and that the 
escape fraction for ionising radiation from the halos is unity. To mimic the strong radiative 
feedback in H II regions, they implemented a "volume exclusion effect" which disabled 
sources within already ionised regions. The results from their analysis indicated that in 
the standard ACDM cosmology, a sufficient number of population III sources turn on to 
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fuUy reionise the volume by 2; ~ 18. By merging the subsequent ionisation history from 
population II sources (adopted from Paper I), they were able to compute a total optical 
depth to Thomson scattering of Te ~ 0.14, in reasonable agreement with the WMAP result. 

Encouraged by the results from Yoshida et al. (2003b), the aim of this paper is to 
use the numerical approach outlined in Paper I to perform a more rigorous study of the 
parameter space associated with the ability of population III stars to reionise the Universe. 
In particular, by tracking individual clouds of dense molecular gas directly, we are able to 
implement simulated feedback effects explicitly. We explore a number of plausible models in 
which we incorporate feedback effects in an ad hoc fashion and examine the corresponding 
ionisation history of the Universe. By coupling this history with the results from the 
population II analysis conducted in Paper 1, we are able to assess the overall success of each 
model in terms of the WMAP and SDSS observations. 

This paper is organised in the following manner. In §2 we describe the cosmological 
simulation used in our analysis. Then in §3 we describe our methodology for source selection 
and present a detailed discussion of how we approximate related feedback effects. The 
results of the simulation are presented in §4. Included in this section is a discussion related 
to the reionisation history of the Universe from the combined effect of population II and 
population III stars. Finally in §5, we provide a summary and state our conclusions. 

2. UNDERLYING COSMOLOGICAL SIMULATION 

The underlying simulation used in this paper was performed using the parallel 
Tree-PM/SPH solver GADGET2, in its fully conservative entropy form (see Springel & 
Hernquist 2002). The simulation follows the non-equilibrium reactions of nine chemical 

species (e~, H, H"*", He, He"*", He"'"+, H2, H2 , H~) using the reaction coefficients compiled 
by Abel et al. (1997). The cooling rate of Galli & Palla (1998) is employed for molecular 
hydrogen. In this paper, we will focus on a ACDM model with matter density Vim = 0.3, 
baryon density Qj, = 0.04, cosmological constant = 0.7 and expansion rate at the present 
time Ho = 70 km s~^ Mpc^^. The spectral index of the primordial power spectrum has 
been set to = 1 and we normalise the fluctuation amplitude by setting cxg = 0.9. Further 
details of the simulation and methods for generating initial conditions are given in Yoshida 
et al. (2003c) and Yoshida et al. (2003d), respectively. 

The particular simulation we use employs 2 x 324^ particles in a cosmological volume 
of 1 Mpc on a side, corresponding to a mass resolution of 160 Mq and 1040 Mq for gas and 
dark matter particles, respectively. Convergence tests using higher resolution simulations 
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indicate that the mass resolution adopted here is sufficient to follow the coohng and collapse 
of primordial gas within low mass (~ IO^Mq) halos. 

Owing to the relatively small size of our simulation box, the fundamental fluctuation 
mode starts to become non-linear at 2; ~ 16, forcing us to stop our analysis at this redshift. 
We discuss the hmitations of such a small volume in the context of predicting cosmic 
reionisation in §4. 

3. SOURCE DEFINITION AND FEEDBACK EFFECTS 

3.1. Source Definition 

Our approach for defining sources relies on identifying virialised halos which contain a 
reservoir of cold and dense gas that undergoes star formation. This requires first locating 
dark matter halos via a friends-of-friends (FOF) algorithm with a linking parameter 
b = 0.164 in units of mean particle separation, and discarding groups with fewer than 100 
dark matter particles. For each FOF group, we compute a virial radius i?vir defined as the 
radius of the sphere centered on the most bound particle of the group having an overdensity 
of 180 relative to the critical density. Within this radius we then measure the mass of gas 
which is cold (T < 500 K) and dense (rin > 500 cm~^). The clouds of molecular gas in these 
reservoirs is expected to rapidly cool and exceed the characteristic Jeans mass Mj ~ 3000 
M0 for the typical temperature T 200 K and density n-n ~ 10^ cm^^ of the condensed 
gas. The halos which contain these reservoirs of cold gas are therefore considered to be sites 
of active star formation and serve as potential locations for population III sources. In the 
top panel of Figure 1 we show how cold gas accumulates in star-forming halos as a function 
of redshift. Here, star-forming halos {solid circles) are plotted at each simulation snapshot 
(separated in time by 3 x 10^ years) and serve to illustrate how cold gas accumulates 
within individual halos over time. In the bottom panel of the same figure we show the 
total number of star forming halos as function of redshift in our simulation volume. It is 
important to point out that in carrying out our analysis, we construct a halo merger history 
and use it to carefully track the star-forming gas that is exchanged between the progenitors 
and descendants of each halo. 

For each potential source site, we compute the amount of gas that actually forms stars 
by multiplying the total amount of cold gas in the halo by the free parameter E^f- Here 

can be thought of as a universal star formation efficiency of dense-cold gas as defined 
above to form sources and should not be confused with similar definitions of star formation 
efficiencies based on the total amount of gas in the halo (see, e.g. Cen 2003a,b). In this 
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study, we consider two different values for this parameter: E^f — 0.02, and Egp — 0.005. 
For a given value of ^'sf, we apply a further restriction on the allowable mass of stars 
which can actually form. More specifically, we assume that the primordial IMF is 
moderately top-heavy and yields stars with masses in the range 100-300 M0. We therefore 
first discard all halos with star masses below 100 M0. Above this limit, we allow a range 
of star masses to be incorporated into the halo up to a maximum of 900 Mq under the 
assumption that only a few stars can form at coincidentally close times and thus survive 
the intense dissociating radiation from the first star in the halo (see Omukai & Nishi 1999; 
Abel, Bryan & Norman 2002) and its subsequent explosion (e.g. Bromm et al. 2003). 
Our assumption may therefore be considered as a slightly more relaxed version of the 
"one-star-per-halo" prescription employed in related studies of this topic (see, e.g.. Oh, 
Nollett, Madau, Wasserburg 2001). Nevertheless, it should be noted that the associated 
parameter space which will be explored in this study leads to relatively few multi-source 
halos (M* > 300 M0) under the assumption of a moderately top-heavy IMF and therefore 
leads to similar results. 

Because of their high mass, population III stars are dominated by radiation pressure 
and have luminosities close to the Eddington limit LEdd- This leads to a constant source 
lifetime estimated to be t ~ 0.007Mc^/LEdd ~ 3 x 10^ yr. We use this value as the universal 
lifetime for all our stars and also adopt it as the time step over which we evolve our 
simulation. To compute ionising intensities we first note the results presented in Bromm, 
Kudritzki, & Loeb (2001) who demonstrated that population III stars with masses > 300 
Mq have a generic mass-scaled spectral form that is almost independent of mass resulting 
in H I and He II ionising photon production rates 1.6 x 10^^ s"^ Mg^ and 3.8 x 10^^ s~^ 
Mq-"^, respectively. For 100 Mq stars, the number of ionising photons per solar mass is 
reduced by a factor of ~ 2 for H I and ~ 4 for He II (Bromm ct al. 2001; Tumlinson & 
Shull 2000). Given our assumption of the 100 - 300 Mq mass range for individual stars, 
we compute ionisation rates for halos by interpolating between the values quoted above for 
the 100 Mq and > 300 Mq cases^. In the relatively few halos where the total source mass 
exceeds 300 Mq, we adopt the median values of 1.2 x 10^^ s"^ M"^ for H I and 2.4 x 10^^ 
s~^ Mq^ for He II under the assumption that these halos contain multiple stars with masses 
in the range 100 - 300 Mq. 

The quantity of ionising fiux that can actually escape a halo and enter into the IGM is 



^Note that in our simulations we carry out radiative transfer calculations to track H II and He III regions 

only, with the implicit assumption that for massive stars, the boundaries of the ionisation zones for He II 
and H II coincide. This serves to significantly speed up our calculations by allowing the tracking of both 
zones via a single ray tracing calculation in H II; see §2.1 of Paper I for further details. 
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parameterised by the escape fraction /esc- In this context, /esc is defined to be the fraction 
of ionising photons intrinsic to each source that escape the virial radius of the halo and 
participate in the radiative transfer calculations. It must be noted here that this parameter 
inevitably also carries with it any uncertainties associated with the shape of the galaxy, the 
gas and stellar density profiles, and the IMF and star formation efficiency. At low redshifts, 
observations deduced from z < 3 starburst galaxies (see Heckman et al 2001; Hurwitz et 
al. 1997; Leitherer 1995) suggest small values for /esc around /esc ^ 0.10. However, it is 
important to note that these observational results may be underestimating the true values 
owing to undetected absorption from interstellar components. In fact, observations of 29 
Lyman break galaxies at z ~ 3.4 by Steidel, Pettini, Adelberger (2001) appear to suggest 
much larger escape fractions. Furthermore, it is unclear how the relevant factors associated 
with the escape fraction vary with rcdshift. Theoretical work on f^sc at high redshifts (Wood 
& Loeb 2000, Ricotti & Shull 2000) finds small values, although these estimates are highly 
uncertain owing to the lack of information regarding the IMF and star formation efficiencies 
at early times. In addition, the effect of dust, complex gas inhomogeneity and gas dynamics, 
all of which may strongly influence escape fractions, are not included in these theoretical 
studies. More recently, Whalen et al. (2003) have studied the ionisation environment of the 
first luminous objects using detailed one-dimensional radiation hydrodynamical simulations 
and find that the transition from D to R-type fronts occurs within only ~ 100, OOOyrs 
after the star reaches the zero age main sequence. Because these I-fronts exit the halo on 
timescales much shorter that the stars' main sequence lifetime, their host halos can have 
UV escape fractions of ^ 0.95. 

Given these recent results we employ relatively large escape fractions, /esc = 1 and 
/esc = 0.3, in all the models we consider. The escape fraction of unity is motivated by 
the extreme cases from Whalcn ct al. (2003) while the escape fraction of 0.3 is meant to 
represent an alternate case where the surrounding environment is much more complex and 
restrictive to ionising break-throughs. 

3.2. Feedback effects 

Feedback effects from the first stars inevitably exert prompt and significant impact 
on subsequent star formation. Locally, population 111 stars arc expected to significantly 
affect their environments when they explode as supernovae (e.g. Bromm et al. 2003). This 
"internal" feedback will perturb the gas and effectively regulate subsequent star formation 
in the pregalactic cloud. Feedback from the first stars can also affect structure formation 
on a global scale. This is because the UV photons below 13.6 eV produced by these stars 
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can easily propagate through the Universe and build up a soft UV background that can 
alter the chemistry of distant regions. More specifically, studies have suggested that even a 
feeble UV background can photo-dissociate the fragile H2 molecules necessary for cooling of 
the primordial gas (see, e.g., Haiman, Rees, & Loeb 1997; Haiman, Abel, Rees 2000; Ciardi, 
Ferrara, & Abel 2000; Machacek, Bryan, & Abel 2001). Conversely, it has been argued 
that feedback may actually enhance H2 formation ahead of of the ionisation front (Ricotti, 
Gnedin, & ShuU 2001) or through the boosting of the electron fraction in dense regions 
by X-rays emitted from, for example, distant supernovae remnants (Haiman, Rees & Loeb 
1996; Cen 2003a; but see also Machacek, Bryan, & Abel 2003 who find that the net effect 
of X-rays is rather mild) . 

Presently, the technical difficulties involved with the incorporation of radiative transfer 
calculations directly into cosmological simulations prevent us from studying these effects 
numerically in a fully consistent manner. We therefore adopt a more ad hoc approach 
which attempts to include feedback effects through the incorporation of simple, physically 
motivated criteria which dynamically constrain the amount of stellar mass that forms. In 
the following two sections we will describe our prescription for parameterising internal and 
external feedback effects within the context of our simulation. 



Stars with masses in the range 140 Mq ^ M ^ 260 Mq are expected to die as 
supernovae which release up to 10^^ ergs of kinetic energy into their surroundings (e.g. 
Heger & Woosley 2002). Since the gravitational binding energy for gas in sub-galactic 
halos is Eb ~ 10^^[M/(5 x 10^)]^/^[(1 + z)/20], these explosions can severely disrupt their 
immediate surroundings, possibly driving away a large portion of the gas from the halo (see, 
e.g., Bromm, Yoshida, & Hernquist 2003). To mimic this effect, we introduce an "internal 
disruption factor" Df which will parameterise the amount of gas that gets blown out. More 
specifically, for every halo which turns on with star mass M*, we flag a corresponding 
amount of star forming gas mass Msf (cold gas) equivalent to Df x within that halo as 
unusable for future star formation. The prescription for calculating a new star mass for a 
given halo is then quantified at each time step i in terms of E'sf by the expression. 



where Mgp^ represents the amount of usable star-forming gas mass from the previous step 
and AMgp is the amount of new star-forming gas accumulating in the halo at the current 
step. This prescription requires us to explicitly track the history of gas in each halo directly. 



3.2.1. Simulating internal feedback effects 
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accounting for the accumulation of new star-forming gas by continuing accretion from 
surrounding gas, and also from mergers with other star-forming halos. In this analysis, 
we explore where the disruption factor is small, -D/ = 40, and also where 

the disruption factor is large, Df = 120. We admit that these choices for the disruption 
factor are somewhat arbitrary. However, our aim here is to use Df in conjunction with Esf 
to create a range of different models which span the plausible parameter space associated 
with internal feedback effects. In particular, our choice oi Df = 40 will describe a scenario 
where clumps of condensed gas in the halo only inefficiently absorb the energy released 
from a local supernova (possibly owing to a spatially clustered distribution with a small 
cross section and large column density to the explosion site) and are therefore marginally 
disrupted (i?sF ^ Df = 0.2 — 0.8). Such a low disruption factor may also be motivated 
by the potential positive feedback which can arise from propagating compression shocks 
induced by the exploding star which may accelerate star formation within the halo (Mackey, 
Bromm, & Hernquist 2003; Salvatcrra, Ferrara & Schneider 2003). On the other hand, our 
choice of Df = 120 will simulate the case where energy absorption is efficient and effectively 
disperses a large amount {Esf x Df = 0.6 — 2.4) of the gas. Here, values for Esf x Df 
close to unity correspond to the total disruption of star-forming gas within the halo, while 
values larger than unity are meant to represent an extreme disruption of the gas extending 
beyond the halo. Formally, the latter scenario results in negative star masses for some of 
the halos immediately after a star-forming episode. In these cases we simply set the star 
mass to zero, but continue to track the amount of negative star forming gas mass in the 
halo, which decreases in time as new gas is accumulated. Negative star forming gas, in 
this context, represents a ramification of strong internal feedback effects which not only 
disrupt the available star forming gas in the halo but also inhibit the future accretion of 
neighbouring gas surrounding the halo. In all cases, it is important to keep in mind that our 
parameterisation of the disruption factor Df also carries with it the uncertainty associated 
with whether or not the stars formed actually end their lives as supernovae. 

In Figure 2, we plot all the possible star-forming halos in our simulation volume for the 
four different combinations of Df and Esf we tried. Here, the two dotted horizontal lines 
delineate the range of allowable star masses in halos {solid circles) as described earlier. In 
particular, halos with source masses falling short of 100 Mq {crosses) are excluded while 
those above the 900 M0 limit {open- circles) are included, but with star masses artificially 
set to the upper limit of 900 M0. A comparison between this figure and the top panel in 
Figure 1 clearly demonstrates the effect of including a disruption mechanism in the halos. 
Namely, once halos acquire enough star mass to turn on, the surrounding star forming 
mass in that halo is disrupted to varying degrees and does not monotonically increase as 
a function of time. Consequently, a halo which has turned on at a given time may not 
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have enough undisrupted star-forming gas at subsequent time steps to turn on again. This 
effect is especially noticeable in the right-hand column where the disruption factor is high 
with the most dramatic case indicated in the top-right panel. Here, the combination with 
a high efficiency factor renders most halos incapable of turning on more than once. In the 
left-hand column, the disruption factor is 3 times smaller and most halos are capable of 
turning on repeatedly. This is especially true in the bottom-left panel where the efficiency 
is also small and individual halos can be seen switching on continuously once they surpass 
the minimum star mass threshold. Interestingly, in all cases the most massive halo in our 
simulation volume appears to contain enough star-forming mass at redshifts below 20 that 
the disruption caused by a star mass of 900 switching on at each step does not inhibit 
subsequent star formation. This seems like a plausible scenario for halos with very deep 
potential wells that are especially conducive to star formation. It must be noted that at 
some level there exists a degeneracy between the effects produced by the parameters Df 
and E^F- However, as Figure 2 demonstrates, there arc distinct differences between the two 
parameters in terms of simulating the physical consequences of star formation, especially in 
extreme cases. We therefore continue our analysis with this distinction. 

3.2.2. Simulating external feedback effects 

We now turn to the issue of simulating external feedback effects which have a more 
global impact. The first type of external feedback we consider is related to the "sphere of 
influence" associated with the ionisation zone surrounding a source. The gas within this 
region is expected to be substantially photo- heated to temperatures in excess of 2 x 10^ K 
when ionised by radiation with a spectrum typical of metal-free stars (see Hui & Haiman 
2003). For example, a 7 x 10^ M© halo at z ~ 20 has a binding energy of roughly ~ 10^^ 
ergs, which is considerably less that the total ionising energy produced by a single massive 
star of order ~ 1 x IQ^^^Mstar/'^^^ Mq) erg. As a result, star-forming gas in shallow 
potential wells located within these ionised regions photo-evaporates and fails to produce 
stars. The magnitude of the effect is difficult to gauge without a self consistent treatment 
of radiative photo-heating effects linked directly to the hydrodynamics of the gas. Since 
the radiative transfer and gas dynamics are not coupled in our calculations, we can only 
simulate the effect in an ad hoc manner. In particular, we would like to qualitatively mimic 
the inhibiting effects that photo-heating may have on star formation in halos as function 
of the depth of their potential wells, or masses, as seen in the simulations of e.g. Bromm 
et al. (2003) and Whalen et al. (2003). For this, we implement an "ionisation exclusion 
effect" which will act to inhibit sources from turning on before a certain amount of time 
has passed since their gas was ionised. We express the recuperation time trecup in terms of 
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the free fall time associated with each source cell. In a spherical potential, the free fall time 
can be expressed as, 



where G is the gravitational constant and p is the density. We note that free fall times 
are calculated only for the overdense cells which actually contain sources, otherwise this 
time scale is not meaningful. As a first order approximation, we use this expression to 
compute free fall times in source cells according to their matter densities. We can then 
express the recuperation time in terms of the free fall time trccup = -^^rccup^fr- Here the free 
parameter K^ccnp can be thought of as the relative fraction of a free fall time that must 
pass before a source halo can recuperate from the negative effects of photo-evaporation. 
This method requires us to carefully track when source cells become ionised (either by their 
own radiation or by the radiation of neighbouring sources) and subsequently increment 
the amount of time that passes in these cells at each time step. It should be noted that 
we also track the migration of sources into different source cells during the course of a 
particular exclusion interval so as to properly follow the total wait time associated with 
each source halo. Owing to the large overdensities in source cells, recombination times are 
short and, as a result, it is possible for a given source cell to undergo multiple ionisations 
and inhibited phases. On a general level, this effect is very similar to the "volume exclusion 
effect" adopted by Yoshida et al. (2003b) although our method has the convenient feature 
of naturally gauging the resiliency of each halo to be able to recuperate from photo-heating. 
To explore both a weak and strong case of this effect, we consider recuperation parameters 

-^recup = 1/3 and -Krecup = 2/3. 

Next, we turn to photo-dissociation of molecular hydrogen by the LW radiation from 
the first stars. While photons with energies above 13.6 eV are likely to be completely 
absorbed by the gas surrounding the sources, those with energies in the LW bands can 
easily travel unimpeded through the neutral IGM and build up a uniform UV background 
radiation field (Dekel & Rees 1987; Haiman, et al. 2000; Omukai & Nishi 1999). We 
emphasise that the net effect of this radiation is not to completely inhibit primordial gas 
cooling, but it raises the minimum mass scale of the halos in which the gas can cool (Haiman 
et al. 2000; Machacek et al. 2001). Recently, Yoshida et al. (2003c) studied the effect of 
LW backgrounds on the formation of primordial gas clouds and found that owing to the 
photo-dissociation of H2, gas cooling is suppressed for radiation with intensity J21 > 0.01 (in 
units of 10^^ erg s~^ cm~^ Hz~^ str~^). They also implemented a technique to compute H2 
column densities around virialised regions making it possible to estimate gas self-shielding 
factors in their maximal limits. With the help of self-shielding, primordial gas in large halos 
can cool efficiently even when exposed to background radiation with intensity J21 = 0.01. 




(2) 
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While the two extreme cases studied by Yoshida et al. (2003c), an optically thin hmit and 
maximal gas self-shielding, should bracket the true effect of the LW background, a more 
careful study is required to obtain an accurate estimate of the impact of this radiation. 

Given the uncertainties, we continue with our approach of adopting a simple 
parameterisation aimed at capturing the qualitative aspects of the effect. In particular, 
we first note that LW radiation is only capable of affecting the abundance of molecular 
hydrogen in gas and does not significantly affect the thermal state of the gas. Assuming 
equilibrium, the amount of molecular hydrogen in a gas cloud is inversely proportional to 
the photo-dissociation reaction coefficient given by. 



where Jlw is the radiation intensity at hu = 12.87 eV and -ft'shicW is the effective shielding 
factor. Neglecting the evolution of the gas density and temperature, and noting that the 
cooling rate scales just as the number of hydrogen molecules, the amount of gas that can 
cool should be proportional to l/( ^LW-f^sMeid)- Given this linear scaling we introduce a 
modulating function. 



from which the amount of cold gas available for star mass production is computed after the 
effects of photo-dissociation are taken into account. To determine the fitting parameters 
a, 6, and Xo, we refer to Yoshida et al. (2003c) who point out that LW radiation is 
unimportant for Jlw < 10"^'^, whereas for Jlw > 10~^^, the gas cannot cool by molecular 
hydrogen cooling in nearly all the halos in the optically thin limit (-^'shielding = !)• Thus by 
setting a = 10^^, b = 1, Xo = 10^^^, and restricting the maximum value of Flw to 1, we are 
able to mimic this behaviour. The fact that the relevant range of intensities is small lends 
support for our use of a linear approximation. 

Assuming that an effective screen owing to abundant neutral hydrogen blocks photons 
in the Lyman-series lines from all sources at redshifts above ^max (Haiman et al. 1997), 
we can approximate Jlw at each redshift Zi by summing the emissivity in the LW band 
from all previously active sources up to z^g^. For our analysis, we set z^^^.^ = Zi^i and 
approximate Jhw{zi) as arising from the unattenuated radiation field purely from the set 
of sources in the previous time step. Our assumption should be reasonable given the large 
amount of neutral hydrogen present around the time when Jlw becomes appreciable and 
the fact that photons with E — 13.6 eV will be redshifted out of the relevant energy range 
(11. 18-13. 6eV) after ~ 18% change in redshift. Since the luminosity in the LW band per 
unit stellar mass for very massive stars is Llw — 3 x 10^^ erg s~^ Hz~^ M©^) with only a 



^diss — 1-38 X 10 JLW-f^shield, 



(3) 



Flw — a( JLW-f^shield — ^o) + b, 



(4) 
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weak dependence on the stellar mass, our expression for Jhwi^i) becomes, 

JLwizi) = ^ Llw X Mf*(^i_i) X Ati X V^-i{zi), (5) 

where c is the speed of hght, M^°*{zi-i) is the total star mass which was switched on at 
redshift step Zi-i, Ati is our time step of 3 x 10^ yr, and V^^{zi) is the proper volume of 
our simulation box at Zi. 

Owing to the highly uncertain nature of shielding processes, we refrain from developing 
a complex parameterisation of this effect (see, however, §7.2 of Yoshida et al. 2003c). 
Rather, we adopt the conservative value of i^^ghieid — 0.10 as our maximum shielding 
effectiveness for halos with M^f > 10^ Mq and linearly scale this value with mass such that 
^shield = 1 for MsF < 10^ Mq (corresponding to zero shielding). We apply this level of 
shielding in all our models. 

Finally, we turn to the effect of metal pollution from the hfe-cycles of the first stars. 
Metals can cool gas inside halos to temperatures lower than can be achieved by H2 
cooling, making possible the formation of less massive ordinary population II stars. The 
critical transition between the two populations should occur after a certain level of metal 
enrichment has occurred. Assuming that the metallicity in the IGM is directly linked to the 
amount of gas that has been incorporated into population III stars, one can parameterise 
when this transition should occur in terms of the fraction of total gas that forms massive 
stars. Oh et al. (2001) estimate that the critical transition should occur when a fraction 
of 3 X 10~^ — 1.2 X 10~^ is formed into massive stars in the range 150-250 M0, which are 
believed to be the main contributors of metals. It should be noted, however, that this 
estimate assumes that metals produced by massive stars are uniformly mixed throughout 
the IGM, which is unlikely. If this is not the case, and metals remain preferentially near 
sites of halo formation, then we would expect that the critical fraction is potentially much 
lower than the value reported by Oh et al. (2001). In addition, stars with masses larger than 
250 Mq were excluded under the assumption that their contribution to metal enrichment 
should be small if they end their lives as black holes. However, the amount of mass these 
stars shed before they collapse into black holes is still relatively uncertain and may also 
reduce the fraction quoted above. 

Given the large uncertainty associated with this fraction, we allow all our models to 
continually form population III stars all the way to 2; = 16 where we stop our simulations 
owing to the growing non-linearity of the fundamental fluctuation mode. For all the models 
which we consider, the fraction of total gas mass that is converted to stars Fconv by 2; = 16 
falls short of the lower bound of the critical range quoted by Oh et al. (2001). In an effort 
to study more protracted population III epochs, we will employ simple extrapolations of 
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the model results to lower redshifts. This will allow us to assess the impact of population 
III stars for larger Fconv values more consistent with the critical range quoted in Oh et al. 
(2001) (see §4.3). 

3.3. Models 

Our goal of identifying the important features which are associated with the ability of 
population III stars to reionise the Universe leaves us with a rather large parameter space 
to explore. Namely, even neglecting uncertainties associated with the LW background and 
metal enrichment, we are still left with 4 free parameters: /esc, -E'Sf, -Df, and -ft'recup- Our 
approach of considering both a weak and strong case for each of these parameters thus leads 
to 16 separate models which we list in Table 1. 

In addition, we include a further model M9 which will represent our most restrictive 
case where we allow only a single 200 star to form in each halo which has at least 4 x 10^ 
Mq of star forming mass {E^f = 0.005). This model is meant to be more consistent with 
the picture described by Abel et al. (2002) who have used high resolution hydrodynamical 
simulations to show that a pre-galactic halo tends to form only a single collapsing core in 
its center without renewed fragmentation. In an effort to simulate the ultimate inhibition 
of subsequent star formation in these halos, we further restrict the model to allow only one 
episode of population III star-formation per halo without the possibility of forming new 
metal-free in the same halo or its descendents at a later time. In this case, once a single 
source turns on within the halo or its gas becomes ionised by a neighbouring source, we 
inhibit it from hosting any future metal-free stars {Df = oo, -ft'rccup = oo). It is important 
to point out that the total number of ionising photons released by a single 200 M0 star in 
the course of 3 x 10^ yrs easily exceeds the amount of baryons in even our most massive 
halos and therefore justifies our use of relatively large escape fractions in single-star halos 
(see Whalen et al. 2003). 

The main goal of our analysis will be to explore the effectiveness of population III 
stars as ionising sources in the early Universe. By considering a plausible set of models 
which incorporate feedback effects to varying degrees, we develop a better qualitative 
understanding of the complex interplay between these processes and reionisation. Our 
inclusion of an ultra-restrictive model will allow us to loosely constrain the minimum 
effectiveness of population III reionisation. 
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4. RESULTS AND DISCUSSION 

The reionisation history of H I and He II is followed in each of our models using the 
radiative transfer code described in Paper I. Briefly, radiative transfer is performed on 
a Cartesian grid with 200'^ cells using an adaptive ray-casting scheme (Abel & Wandelt 
2002). Density fields, clumping factors, and source characteristics are taken from outputs 
of the hydrodynamic simulations; thus we neglect the dynamical feedback of the radiation 
on structure formation^. However, to account for the neglected photo-heating of the IGM, 
we compute recombinations in our radiative transfer calculations according to an artificially 
raised temperature of T = 1.5 x IC^ K. This temperature corresponds to a reasonable 
estimate of the IGM temperature after reionisation by sources with spectra typical of 
metal-free sources (Hui & Haiman 2003). If the photoionised gas has a higher temperature, 
recombinations occur at a slightly slower rate (the recombination time tree oc T°'^) and 
ionisation zones grow marginally faster. For lower temperatures, the opposite is true. 

In the context of reionisation studies, cosmic variance limits our ability to capture both 
a representative mass function for source halos and also a reliable picture of the distribution 
of gas in the Universe. Since our simulation box has a comoving length of only 1 Mpc, we 
will not be able to make very strong statements regarding the precise moment of cosmic 
reionisation. However, as stated earlier, our goal is only to assess the relative properties 
of a set of physically motivated models for the impact population HI sources can have on 
reionisation. Having stated this, we note that the same small box simulation was analysed 
in Yoshida et al. (2003a) and was found to contain halo abundances for M > 7 x 10^ 
Mq that were in reasonable agreement with the Press-Schechter mass function between 
17 < 2; < 25 (see also Jang-Condell & Hernquist 2001). The incomplete sampling of the 
halo mass function due to the finite box size is appreciable only at 2; > 30. However, as 
we shall show, star formation in our models is not dominated by rare massive halos and 
therefore simulating larger volumes should not significantly alter our results. 



4.1. Star Formation History 

In Figure 3, we show the evolution of the star formation rate (SFR) density for 
population HI stars for all our models. Here the SFR is computed on the fly as the 
simulation regulates the amount of star mass that actually turns on after incorporating the 



^See Sokasian, Abel & Hernquist (2001) for a more detailed description of how we integrate radiative 
transfer calculations with existing outputs from cosmological simulations. 
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effects of feedback. By z — 16, all but our most restrictive model M9 appear to reach a 
SFR near 10~^ yr~^ Mpc~^. The similarity in the SFRs between the different models 
brings attention to the self-regulating consequences of feedback. For example, when the 
disruption factor is large, less star formation takes place locally, which in turn minimises 
negative feedback in the form of dissociating radiation and photo-evaporation. As a result, 
star formation is promoted externally and the global SFR is regulated. In a similar fashion, 
models with a low escape fraction limit the effect of photo-evaporation externally and 
therefore are able to have a similar level of star formation as models with a large escape 
fraction. In §4.3 we shall take a more detailed look at dynamic feedback and attempt to 
quantify the effect. 

Nevertheless, it is important to point out that a SFR of 10~^ Mq yr~^ Mpc~^ is 
roughly an order of magnitude lower than the "normal-mode" (population II) star formation 
predictions from Springel & Hernquist (2003) and Hernquist & Springel (2003) at the 
same redshift. Here, "normal mode" refers to star formation occurring in larger mass 
systems where gas cooling takes place via atomic hydrogen and helium transitions and the 
only form of regulating feedback is supernovae. Besides the difference in amplitude, the 
SFR in our population 111 models appears to have relatively "flatter" shapes in the range 
16 < ;2 < 20 compared to the "normal- mode" SFR of Hernquist & Springel (2003). This is 
not surprising given the stronger regulating processes associated with population HI star 
formation. The overall shape and the amplitude of the SFR are in reasonable agreement 
with the semi-analj^ic prediction of Yoshida et al. (2003c, their Fig. 18). It is important 
to note that population HI star formation at z ^ 20 may have a significant impact on 
subsequent "normal-mode" star formation occurring in small proto-galaxies which arc 
susceptible to the negative feedback effects of photoheating. Because these effects were not 
included, it is not surprising that the population 111 SFR presented here does not converge 
with the "normal- mode" SFR of Hernquist & Springel (2003). Nevertheless, if the epoch 
of population 111 stars ends abruptly, as we have assumed here, then the transition itself 
should also be abrupt. A very interesting issue for future work, therefore, would be to 
address how the transition from population HI stars to ordinary star formation occurs. The 
discontinuity presented here represents an additional uncertainty of our analysis only when 
we attempt to compute a complete history of reionisation which spans both population II 
and population HI stars. 
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4.2. Global lonisation Fractions 

The evolution of the global ionisation fraction in the simulation volume is a useful 
tool for characterising differences between the various models. In Figures 4 and 5 we show 
the evolution of the volume-weighted ionisation fraction down to 2; = 16 for H II and He 
III, respectively. In each panel, we plot both the /esc = 1-0 {solid-line) and /esc = 0.3 
{dashed-line) cases for the corresponding model. Here, we note that with the exception of 
our most restrictive model M9, all the /esc = 1-0 models are able to significantly ionise 
their volumes by ^ = 16. Interestingly, some of the models with higher emissivity are also 
able to substantially doubly ionise helium although the corresponding fractions are much 
less than in H II. This is true despite the fact that the ionising spectrum of a typical 300 
M0 metal-free star produces roughly ~ 3 times more ionising photons per helium atom 
compared to hydrogen for cosmic abundances. The difference in the resultant ionisation 
fractions owes to the fact that the recombination rate of He III is much larger relative to H 

II. This large recombination rate is also responsible for making He III ionisation evolution 
relatively less smooth, because relic He III regions are much quicker to recombine. 

In all our models, the initial rise of the ionisation fraction appears to be closely 
correlated with the exponential growth of star-forming halos in the volume. Interestingly, 
however, our highest emissivity models (Mlb, M2b, M3b) appear to show some degree 
of divergence from this rate at later times. In particular, the rise of the ionisation 
fraction becomes shallower once the ionisation fraction exceeds ~ 60%. The mostly likely 
explanation for this behaviour is that dynamic external feedback becomes particularly 
inhibitive once a substantial portion of the volume has become ionised and many sources 
become involved in the process. We will explore this effect in more detail in the next section 
where we try to assess the impact from this dynamic component of the feedback in terms of 
the amount of star mass that actually turns on at each time step. 

A summary of the ionisation results at 2; = 16 is given in Table 2 which lists for each 
model: the fraction of total gas mass converted to stars -Fconv, the cumulative number of 
ionising photons released per atom, and the volume ionisation fraction for H II and He 

III, respectively. Here we point out that while our more restrictive models do not achieve 
as high of an ionisation fraction as the models with larger emissivities by z — 16, they 
convert less gas mass into stars and therefore should have polluted their environment 
relatively less with metals. As a result, we would expect that metal-free star formation 
should persist longer in these models. It is particularly interesting to compare the results 
for our most restrictive model M9 to our most liberal model Ml. We note in particular 
that the realisation of a "one-star-per-halo" for metal-free stars as defined in our M9 model 
has such a slow rate of growth for its ionising emissivity that recombinations will likely 
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prevent the ionisation fraction from increasing significantly even if we had extended this 
model down to lower redshifts. However, wc point out that this model is relatively more 
efficient in terms of the resultant ionisation fraction per ionising photon than the Ml 
model with larger emissivity. This feature can be attributed to the fact that models with 
higher emissivity suffer cumulatively more recombinations as a result of ionising substantial 
portions of their volumes at very early redshifts and keeping them ionised down to lower 
redshifts. It should be noted, however, that this difference will be shghtly diminished by 
the fact that ionisation fronts in the more restrictive M9 model have difficulty breaking 
through the dense environments surrounding the sources. This causes the M9 model to 
suffer relatively more recombinations per unit volume than models with higher emissivity 
whose ionisation zones expand into the less dense and more voluminous portions of the 
IGM. This morphological effect also helps to explain why models M4-M9, which have not 
yet experienced a slow-down in the rate of ionising photon production (owing to external 
feedback factors), exhibit ionisation fractions for their /esc = 0.3 case that arc systematically 
lower than 0.3 times the ionisation fractions predicted in the corresponding /esc = 1-0 case 
(see Table 2). 

A visual illustration of the reionisation process is shown in Figure 6 where we plot a 
series of projected slices through the simulation volume from the Mia model (other models 
exhibit a similar topological evolution in their ionisation structures). In each panel, a 0.25 
Mpc slice is projected in both density and ionisation fraction. Prom the plots one can 
follow the growth of the ionisation zones {blue) around the first stars as they turn neutral 
gas {yellow) into highly ionised regions (light-blue). The plots clearly show the preferential 
advance of ionisation fronts into the undcrdcnsc regions, responsible for the steep rise seen 
in the corresponding volume-weighted ionisation fractions. 

To understand which halos are contributing the bulk of the ionising radiation, we plot 
in Figure 7 the fraction of ionising flux released as function of halo mass for three redshift 

ranges associated with the Mlf, run (other models had similar results). Here we see that 
through the redshift ranges 21 < z < 23, 19 < 2; < 21, and 16 < z < 19, the bulk of 
ionising radiation is consistently released from halos with masses between ~ 2 x 10^ Mq 
and ~ 6 X 10^ M0. The fact that these "mini-halos" arc capable of substantially ionising 
the IGM appears to contrast with the findings of Ricotti ct al. (2002a) who conducted 
similar 3D cosmological simulations which included dynamically linked radiative transfer 
calculations (see Ricotti et al. 2002b for simulation details). They showed that negative 
feedback prevents the size of H II regions from exceeding the size of the dense filaments; 
hence, "mini-halos" are not able to rcionisc the voids. More specifically, the Ricotti et al. 
(2002a) simulations showed that positive feedback from the accelerated formation rate of 
molecular hydrogen in front of ionisation fronts within these filaments provided an effective 
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shield preventing the destruction of molecular hydrogen in other star forming regions. On 

the contrary, when the ionisation fronts expand beyond the filaments, densities drop and 
the shield disappears allowing negative feedback from LW radiation and photoevaporation 
to dominate. Ricotti et al. (2002a) thus concluded that the volume filling factor of the H II 
regions associated with these "mini-halos" remains small and cannot reionise the universe. 

Although none of the models considered in our analysis are able to fully reionise the 
universe before the onset of population 11 stars, they do ionise a substantial portion of 
the IGM. This discrepancy most likely results from the way in which negative feedback 
is treated in the simulations. More specifically, by dynamically linking approximate 
radiative transfer calculations directly with cosmological simulations, Ricotti et al. (2002a) 
are able to simulate positive feedback processes preceding H II regions which effectively 
reduces negative feedback from the dissociation of molecular hydrogen. However, when 
the ionising emissivity of "mini-halos" is large, ionisation fronts expand outside of the 
filaments and negative feedback from dissociation and photo-heating begins to dominate. 
Since our analysis does not allow us to couple radiative transfer effects directly into the 
cosmological simulations, we are unable to include positive feedback effects such as the one 
mentioned above. Nevertheless, the models we consider are able to easily force ionisation 
fronts beyond the dense filaments and therefore we do not expect positive feedback to be 
important in these cases. Once ionisation fronts expand into the IGM, negative feedback 
effects kick in through our ad-hoc implementation of the photo-evaporation effect and the 
LW dissociating background. As we shall show in the following section, such feedback 
works to regulate further star formation, however, unlike Ricotti et al. (2002a) we do not 
find that it fully suppresses star formation. It is important to note that our simulations 
included an approximate treatment which allowed halos to shield themselves against 
dissociating radiation in proportion to their masses. In contrast, self-shielding was excluded 
in the simulations of Ricotti et al. (2002a) and may have resulted in an overestimation of 
negative feedback, plausibly causing much of our discrepancy with the results by Ricotti 
et al. (2002a). Furthermore, it should be pointed out that the discrepancy may also be 
related to a number of other differences between the two simulations such as the differing 
methodologies used in the calculation of the star forming mass in the halos as well as 
differences in the choice of the IMF. 



4.3. Dynamical Feedback from the LW Background and Photo-evaporation 



In Figures 4 and 5 we showed the evolution of the global volume-weighted ionisation 
fraction for each model. While useful for making comparisons between the net ionising 
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ability of the models, the figures do not reveal the complex role dynamical feedback 
may have played in the evolution. In this section we attempt to study the effectiveness 
and subsequent evolution of this feedback more directly. In particular, we focus on 
tracking the modulating behaviour of our implementation of the LW background and the 
photo-evaporation effect (or in terms of our parameterisation, the "ionisation exclusion 
effect"). More specifically, by tracking the ratio of the star mass that actually turns on at 
each time step M°^^j. to the star mass originally available M°^lf (computed in the absence 
on any external feedback), we can directly study the net impact from both these effects. 
In Figure 8 we show the ratio M°^^^/M°^^ at each redshift step as a histogram for all 
our models. The plots reveal sharp swings mainly owing to the discrete timesteps used. 
However, these swings are enhanced at early times, refiecting the strong nature of our 
"ionisation exclusion effect" when localised to a small region with only a few sources, which 
is the case when the first cluster of sources turns on. As more sources are invoked at lower 
rcdshifts, this ratio begins to behave more smoothly. 

In Table 3 we list for all the models the mean value of M°^/M°^lf weighted by 
the original star mass between the time the first source turns on and z = 16. This 
mass-weighted mean quantifies the relative effectiveness of external dynamical feedback. 
Note, in particular, how models in which we have incorporated strong local feedback effects 
have correspondingly low values for this mean ratio (compare, for example, models Ml-*->M2 
and M3<->M4). This is a clear indication of the self- regulating feature of the population 
III era: when local negative feedback is strong, the global component of negative feedback 
is reduced, allowing for more sources to turn on in distant regions. This feature helps to 
explain how such a large set of model realisations leads to the relatively similar ionisation 
histories seen in Figures 4 and 5. 

Bearing in mind that the net effect of the LW background is to systematically reduce 
Mg°^/Ms°"f in proportion to the star mass that froms, one can use the mean of this ratio 
to isolate and assess the relative impact of the photo-evaporation effect between the models 
M1^M5, M2^M6, M34->M7, and M4^M8. Note in particular, the significant drop in 
this ratio from M1^M5, which gives an indication that external feedback may play an 
important role in scenarios where star formation is efficient and internal feedback is low. 

Interestingly, the effectiveness of external feedback is strong by 2; ~ 16 even in model 
M9 which turns on very few sources. This is most likely owing to the fact that these 
sources reside in only the most clustered regions containing massive enough halos capable 
of hosting sources in this ultra-restrictive model. The close proximity of the sources to one 
another in these regions therefore keeps external feedback in the form of photo-evaporation 
particularly effective in inhibiting a fair portion of the sources from turning on. We mention 
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that, under certain conditions, weak ionising radiation can promote molecular hydrogen 
formation and enhance primordial gas cooling, hence causing possibly positive feedback 
effects (Haiman, Rees, & Loeb 1996; Kitayama et al. 2001). We do not consider the effects 
in our models and thus all the models may be restrictive in this sense. 

4.4. Electron Optical Depth: Population II and III sources 

The recent tentative measurement of polarisation in the CMB by the WMAP sateUite 
(Kogut et al. 2003) indicates that reionisation occurred earlier than implied by the SDSS 
quasars alone (e.g., Spergel et al. 2003). The optical depth to Thomson scattering is 
0.13 < Te < 0.21 corresponding to an instantaneous reionisation epoch between 14 < z < 20. 
In Paper 1, we showed that stellar sources similar to those seen in local galaxies; i.e. 
population II type stars, are unable to reionise the Universe by such high redshifts. 
Meanwhile, the analysis in this paper indicates that metal-free (population III) stars could 
have significantly reionised the Universe at 2; > 16. Cen (2003a,b) and Wyithe & Loeb 
(2003a,b) explored a large number of models of reionisation history using semi-analytic 
methods and argue that, in some specific cases, the Universe could be reionised twice. 

In an effort to study this scenario more quantitatively, we would like to couple the 
ionisation histories from both types of sources and compute the resulting integrated electron 
optical depth. Unfortunately, owing to the high resolution necessary for carrying out 
simulations involving population III stars, we are unable to directly include these sources 
in the larger simulation volumes necessary for studying cosmic reionisation by population 
II sources (galaxies). We are therefore forced to adopt the more approximate approach of 
linking together the ionisation histories from the two sets of simulation volumes. More 
specifically, our approach relies on using the results of our population III simulations to 
estimate the volume-weighted ionisation fraction at redshifts around which it is reasonable 
to assume that star-forming halos have exhausted their ability to form additional metal-free 
stars. At this point in time, we then re-run the /esc = 0.20 population II simulation from 
Paper I (conducted in a 14.3 Mpc comoving box) with the initial conditions that the IGM 
was uniformly ionised to the same ionisation fraction (both in hydrogen and helium) as 
estimated by our population III simulations and also uniformly heated to T = 1.5 x 10^ K. 
Our use of this particular model for population II sources is motivated by our analysis in 
Paper I where we found that given the amplitude and form of the underlying star formation 
predictions, an escape fraction near 20% is most successful in terms of statistically matching 
observational results from the z — 6.28 SDSS quasar of Becker et al. (2001) (see Paper I for 
further details). By linking together the two simulations in this fashion, we are assuming 
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that the onset of population II begins exclusively at redshifts below the transition redshift 
and that there are no population III stars thereafter. Since the relative contribution 
predicted for population II sources in the simulations is small above the transition redshifts 
which we will consider, we can reliably separate the two contributions in this manner. 

Given the fact that most of our population III models ionise the volume to a similar 

extent hj z = 16, we are motivated to simplify our analysis by grouping together these 
models into a single model by averaging all the results at each redshift. We perform this 
averaging for models M1-M8, but continue to keep model M9 distinct as it represents 
our ultra-restrictive case which predicts significantly lower ionisation fractions as a result 
of invoking only a single 200 Mq star once per star forming halo. By z = 16, at which 
point we stop our population III simulations, our combined model M1-M8 with /esc = 0.30 
(/esc = 1) has converted a fraction Fconv = 1-53 x 10~^ (i^conv = 1-26 x 10~^) of its total gas 
into stars. This fraction falls short of the lower limit estimated by Oh et al. (2001) for the 
transition to the population II epoch (~ 3 x 10~^). In an effort to examine a more prolonged 
population 111 epoch culminating with gas-to-star conversion fractions more in-line with 
the estimates from Oh et al. (2001), we apply a simple extrapolation of all the simulation 
results to redshifts below z = 16. More specifically, we extrapolate the simulations using 
the associated rate of change in the models between z = 18 and z = 16. For the combined 
M1-M8 model with /esc = 0.30 (/esc = 1), we extrapolate the population III epoch down 
to z — 14.5 and z — 12.0 resulting in gas conversion fractions: Fconv — 2.20 x 10~^ 
(Fconv = 1.90 X 10-5) and F^onv = 3.32 x 10'^ (Fco„v = 2.86 x 10"^), respectively. In the 
case of model M9 with /esc = 0.30 (/esc = 1), even an extrapolation down to z = 12 only 
yields a conversion fraction of -Fconv = 5.10 x 10"^ (-Fconv = 3.94 x 10"^), well below the 
estimate from Oh et al. (2001). However at z <^ 12, population 11 stars which form via 
atomic line cooling in massive halos begin to dominate the ionising cmissivity relative to the 
population 111 stars incorporated in M9. We therefore adopt z = 12 as our epoch transition 
redshift for this model keeping in mind that the exclusion of population 111 sources below 
z — 12 slightly underestimates the resultant rise in the global ionisation fraction. 

In Figure 9 and 10, we plot the evolution of the H II volume-weighted ionisation 
fraction (top-panel) and the corresponding optical depth to Thomson scattering Tg between 
the present and redshift z {bottom-panet) for the models discussed above with /esc = 0.3 
and /esc = 1.0, respectively. Here the vertical dashed-dotted line represents the redshift 
(z — 16) where we stop our population HI simulations. The meaning of the various curves 
is summarised in the caption to Figure 9. The electron optical depth was calculated from 




(6) 
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where gt — 6.65 x 10~^^ cm^ is the Thomson cross section and ne{z') is the mean electron 
number density at z' . In the calculation of we use mass-weighted ionisation fractions 
and assume that helium is singly ionised to the same degree as the hydrogen component 
for 2; > 3 and doubly ionised everywhere at z < 3. The figure clearly demonstrates how 
population 111 sources may have significantly contributed to the electron optical depth, 
especially in the case where the escape fraction from the mini-halos hosting metal-free stars 
is close to unity. Interestingly, the relatively late onset of population III stars does not cause 
a dramatic recombination era before population II stars begin to continue the reionisation 
process. 

Comparisons with the WMAP results show that only the combined model with 
/esc = 1-0 comes close to the inferred mean value of — 0.17 (dashed-line) with optical 
depths 0.141 and 0.152 for population III extrapolations to 2; = 14.5 and z — 12, 
respectively. In the case of our ultra-conservative M9 model with /esc = 1-0, the resultant 
optical depth is ~ 40% larger than in the case with only population II stars only, but is 
slightly shy of the lower 2a limit of the WMAP measurement. In all cases, however, it 
appears that a large escape fraction is required during the population III epoch in order to 
significantly contribute to the integrated optical depth. 

5. SUMMARY AND CONCLUSIONS 

In this paper we have applied approximate radiative transfer calculations to the 
results of cosmological simulations capable of following early structure formation in an 
effort to study the potential impact that the first generation of metal-free stars may have 
had on the reionisation history of the Universe. By incorporating a series of approximate 
feedback effects, we are able to simulate the complex interrelated processes which may have 
self-regulated subsequent population III star formation. We have examined a set of models 
in which the severity of the effects was varied in order to capture qualitatively how the 
feedback can affect star-formation at early epochs. 

Overall, we find that for a plausible range of values related to our parameterisations, 
population III star formation proceeds in a self-regulated manner. As a result, the net 
impact of these sources in terms of reionisation is fairly insensitive to the relative degree 
of severity of one form of feedback over another. This is true as long internal feedback 
within a halo is not exceedingly large and halos can host multiple generations of population 
III stars. In particular, we find that if we restrict halos of sufficient mass to host only a 
single 200 Mq star once in their existence, then the impact of population III sources on 
reionisation is significantly reduced. 
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In the models where there is no restriction on the number of stars a halo can host and 
the escape fraction is unity, the mean level of ionisation reached by 2; = 16 is 66% and 
29% for H II and He III, respectively. One would therefore expect that at least within the 
context of this description, population 111 stars should have had a significant impact on 
the ionisation history of the Universe. Interestingly, the fact that the helium component 
may also have been ionised to a significant fraction has consequences for the thermal 
history of the IGM (e.g. Hui & Haiman 2003) and could possibly be observable with future 
observations of He II recombination lines in the host halos using the Next Generation Space 
Telescope (Oh, Haiman, & Rees 2001; Tumhnson, Giroux, & ShuU 2001). 

By linking together the simulation results of the population 11 study from Paper 1 
with the results of the population HI analysis presented here, we are able to synthesise a 
full history of how the Universe may have been reionised. We find that even with fairly 
conservative estimates for negative feedback, population HI sources with large escape 
fractions are able to significantly raise the total electron optical depth so that it becomes 
statistically consistent (within the la lower bound) with the WMAP measurement. We 
note that the Thomson optical depths predicted in our analysis appear to be lower than 
those obtained in the semi-analytic treatment conducted by Cen (2003a). Correspondingly, 
population 111 sources in our analysis appear to significantly affect the reionisation history 
of the Universe only at relatively later times {z ^ 19) than the models considered by Cen 
(2003a). Nevertheless, both studies appear to suggest that achieving Tg > 0.17 seems 
quite difficult with reasonable assumptions for the properties of population HI stars. It 
is impractical to make more direct comparisons between the two works given the very 
different approaches employed and the large number of parameters. It should also be noted 
that semi-analytic treatments do no allow for the direct incorporation of negative feedback 
effects and instead absorb these uncertainties directly into star formation efficiency factors. 
As a result, they may oversimplify the true impact of these effects if the amount of star 
mass that turns on in each halo is a sensitive function of halo properties such as proximity 
to other halos, merger history, and total mass (Yoshida et al. 2003c). 

In carrying out our analysis, we were forced to extrapolate a small portion of the 
population 111 ionisation history to lower redshifts owing to boxsize limitations. In doing 
so, we were able to produce gas-to-star conversion fractions which approached the critical 
range described in Oh et al. (2001), estimated to mark the transition to the population 
II epoch. It is important to point out that the simple recipe of Oh et al. (2001) may not 
accurately reflect the transitional stage if metal pollution is not efficient. In particular, in 
the case of an inefficient mixing mechanism, the mctallicity in overdense regions can almost 
reach solar levels even at 2; > 16 for a top heavy IMF, preemptively ending the population 
HI epoch (e.g. Yoshida et al. 2003e). Nevertheless, in all of our models we switched over to 
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a population II epoch at aggressively conservative values below or slightly above the lower 
bound of the estimated range in Oh et al. (2001). Therefore, it is plausible that population 
III star formation continues beyond the transitional redshifts considered here. If that is 
the case, population III stars would continue to significantly ionise the IGM, potentially 
leading to a single reionisation epoch in which population II sources turn on at relatively 
late times and are intense enough to keep the IGM in a highly ionised state. In fact, our 
combined /esc = 1 population III model, when extrapolated down to z — 12, exhibits such a 
single reionisation epoch. Interestingly, the results from this model appear very similar to 
the prediction from our model in Paper I which incorporated an evolving ionising boosting 
factor as an additional multiplicative factor to the intrinsic ionisation rates in population 
II sources. Such an evolution would be indicative of an evolving IMF which became more 
and more top-heavy with increasing redshift. We note however, that such a "vanilla" model 
of reionisation may be inconsistent with the observational constraints from the thermal 
history of the IGM as noted earlier. Future observations of CMB polarisation by the Planck 
surveyor, potentially tractable observations of the 21 cm emission signal from neutral 
hydrogen at high redshifts (see Furlanetto et al. 2003; Ciardi & Madau 2003; Zaldarriaga 
et al. 2003), and searches for Ly-a emission from galaxies at the reionisation epoch (e.g. 
Barton et al. 2003) may help resolve these degeneracies. 

Through our analysis, we have attempted to assess the potential role feedback effects 
may have had on the ability of population III sources to reionise the Universe. While 
our radiative transfer simulations capture nearly all the important physical mechanisms 
operating during reionisation, the precise effects of some processes that are modeled in an 
ad hoc manner remain still uncertain. More accurate studies may soon be forthcoming with 
the incorporation of radiative transfer calculations directly into cosmological simulations 
similar to the one used here. 

Finally, we emphasise that the results here were obtained for a conventional ACDM 

cosmology. If the true power spectrum is actually similar to the running spectral index 
(RSI) model advocated by Spergel et al. (2003), "first star formation" in mini-halos will 
be greatly suppressed (Yoshida et al. (2003a). It is an open question whether the first 
luminous sources at redshifts below 2; ^ 18 in a RSI cosmology can produce a Thomson 
optical depth consistent with the WMAP measurement. 
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TABLE 2 

Summary of model results at z = 16 
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TABLE 3 



Mean fraction for dynamic modulation of star mass 
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Fig. 1. — Top panel: The amount of cold gas (T < 500 K and nn > 500 cm~^) in star forming 
halos as a function of redshift. Star- forming halos are plotted at each simulation snapshot 
(separated in time by 3 x 10^ years). Bottom panel: The total number of star forming halos 
in our simulation volume as a function of redshift. 




Fig. 2. — Total star masses in halos as a function of redshift for different values of the 
disruption factor Df and star formation efficiency E^f (see text for discussion). Horizontal 
dotted lines indicate the range of allowed star masses in halos {solid- circles). Source masses 
falling below 100 Mq mass {crosses) are excluded, while those above 900 Mq {open-circles) 
have star masses artificially set to an upper limit of 900 Mq. 
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Fig. 3. — Evolution of the star formation rate density in each model down to z — 16. 
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Fig. 4. — Evolution of the volume-weighted ionisation fraction for H II in each model down 
to z — 16. 
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Fig. 5.- 



Same as Figure 4, but for He III. 
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Fig. 6. — A series of projected slices through the simulation volume at {top-left to bottom- 
right) z=21.4, 19.8, 18.7, and 17.6. In each panel, a 0.25 Mpc slice (1/4 of the box length) 
from the outputs of the Mia (/esc = 1-0) model is projected in both density and ionisation 
fraction. From the plots one can follow the growth of the ionisation zones (blue) around the 
first stars as they turn neutral gas {yellow) into highly ionised regions {light blue). 
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Fig. 7. — The fraction of total ionising flux released by population III stars in model Ml^ 
as a function of halo mass for redshift ranges 21 < z < 23, 19 < z < 21, and 16 < z < 19. 
Similar results were obtained for models M2-M9. 
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Fig. 8. — The evolution of the ratio of the star mass that actually turns on at each time step 
Mstar ^'^ ^hc star mass originally available Af °taf before the effects of the LW background and 
the "ionisation exclusion effect" were taken into account (see Table 3 for related statistics). 
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Fig. 9. — Top-panel: The evolution of the H II volume- weighted ionisation fraction obtained 
by linking together the ionisation histories of /esc = 0.30 population III simulations with 
/esc = 0.20 population II simulations described in Paper I. Here the vertical dashed-dotted 
line represents the redshift {z = 16) where we were forced to stop our population III 
simulations owing to the small box size. The dotted lines extending below this redshift 
represent extrapolations of the ionisation fraction if the population III epoch had continued 
to lower redshifts (see text for discussion). The models shown are (from top to bottom): 
(1) blue line: the combined (mean) population III simulation results from models Mla-M8a 
(/esc = 0.30) extrapolated down to z = 14.5 before switching to the (/esc = 0.20) population 
II simulation, (2) green line: same as (1) but with the population III epoch extrapolated down 
to z = 12. (3) cyan line: model M9a {ultra-restrictive "one-star-per-halo" model) with a 
population III epoch extrapolated to ^ = 12 before switching to the population II simulation, 
and (4) {red line): the case where only population II sources are present (representing the 
/esc = 0.20 model from Paper I). Bottom-panel: Corresponding optical depths to Thomson 
scattering Tg between the present and redshift z. The inferred value for Te from the WMAP 
measurement is shown {long- dashed line) as is the corresponding la and 2a lower bounds 
{short-dashed lines). 
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Same as Figure 9 but for population III sources with f^sc = 1-0. 



